---
title: "Climate change may alter human physical activity patterns"
subtitle: "Replication code"
author: Nick Obradovich^[Corresponding author.] ^[Belfer Center for Science and International Affairs, Kennedy School of Government, Harvard University.] ^[Media Lab, Massachusetts Institute of Technology.], James H. Fowler^[Departments of Political Science and Medicine, University of California San Diego.]
date: ""
output: 
  html_document:
    fig_caption: true
---

```{r load, echo=T, message=F, warning=F, include=T}
knitr::opts_chunk$set(fig.path='figures/', warning=F, message=F, fig.retina=1, cache=T, autodep=T, echo=F, cache.lazy = FALSE)

rm(list = ls())

library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
#library(kfigr) #figure referencing for markdown
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(gtable) # for plots
library(bit64) # long integers

# historical data and functions
load('./replication.Rdata')
source('./binned_plot_function.R')

# set cores for parallel processing
registerDoMC(20)

# run from command line
# Rscript -e "rmarkdown::render('./main_text_replication_code.Rmd')"
```

```{r figure1, echo=TRUE, fig.width=25, fig.height=9, dev='png', dpi=1200, fig.cap="**Main text Figure 1**"}
# cuttmax regression
tmax.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss, exactDOF=TRUE, psdef=FALSE)
summary(tmax.fe)

heat.fe <- felm(exerany2 ~ cutheat + cutprcp + cutcloud + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss, exactDOF=TRUE, psdef=FALSE)
summary(heat.fe)

# tmax
tmax.p <- ggMarginal(binned.plot(felm.est = tmax.fe,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "Mean Monthly Max. Temp. in \u2103",
             ylabel = "Percentage Point Change in Pr(Physically Active)",
             panel = "a",
             panel.x = -4.5,
             panel.y = -11.5),
           data=brfss,
           x=brfss[, avg.tmax],
           type="histogram",
           margins="x",
           size=20,
           xparams = list(bins = 30, fill = "gold", color = 'gray60', alpha=0.75, size=0.1))

# heat
heat.p <- ggMarginal(binned.plot(felm.est = heat.fe,
             plotvar = "cutheat",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "Mean Monthly Heat Index in \u2103",
             ylabel = "Percentage Point Change in Pr(Physically Active)",
             panel = "b",
             panel.x = -4.5,
             panel.y = -11.5),
           data=brfss,
           x=brfss[, avg.nwsheatindex],
           type="histogram",
           margins="x",
           size=20,
           xparams = list(bins = 30, fill = "gold", color = 'gray60', alpha=0.75, size=0.1))

# prcp
prcp.p <- ggMarginal(binned.plot(felm.est = tmax.fe,
             plotvar = "cutprcp",
             breaks = 2,
             omit = c(0,2),
             ylimit = c(-11.5,1.5),
             xlabel = "Monthly Precipitation Days",
             ylabel = "Percentage Point Change in Pr(Physically Active)",
             panel = "c",
             panel.x = 1,
             panel.y = -11.5),
           data=brfss,
           x=brfss[, prcp.days],
           type="histogram",
           margins="x",
           size=20,
           xparams = list(bins = 30, fill = "gold", color = 'gray60', alpha=0.75, size=0.1))

# humid
humid.p <- ggMarginal(binned.plot(felm.est = tmax.fe,
             plotvar = "cuthumid",
             breaks = 10,
             omit = c(50,60),
             ylimit = c(-11.5,1.5),
             xlabel = "Monthly Average Rel. Humidity",
             ylabel = "Percentage Point Change in Pr(Physically Active)",
             panel = "d",
             panel.x = 15,
             panel.y = -11.5),
           data=brfss,
           x=brfss[, prcp.days],
           type="histogram",
           margins="x",
           size=20,
           xparams = list(bins = 30, fill = "gold", color = 'gray60', alpha=0.75, size=0.1))

grid.arrange(tmax.p, heat.p, prcp.p, humid.p, ncol=4)
```

```{r figure2, echo=TRUE, fig.width=15, fig.height=9, dev='png', dpi=1200, fig.cap="**Main text Figure 2**"}
## young vs. middle vs. old
young.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss[age < 40,], exactDOF=TRUE, psdef=FALSE)
summary(young.fe)

middle.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss[age >= 40 & age < 65,], exactDOF=TRUE, psdef=FALSE)
summary(middle.fe)

old.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss[age >= 65,], exactDOF=TRUE, psdef=FALSE)
summary(old.fe)

## normal vs. overweight vs. obese
normal.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss[bmicat==1,], exactDOF=TRUE, psdef=FALSE)
summary(normal.fe)

ov.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss[bmicat==2,], exactDOF=TRUE, psdef=FALSE)
summary(ov.fe)

ob.fe <- felm(exerany2 ~ cuttmax + cutprcp + cutcloud + cuthumid + cuttrange | mmsa + time + mmsa:season | 0 | mmsa + time
           , data=brfss[bmicat==3,], exactDOF=TRUE, psdef=FALSE)
summary(ob.fe)

## normal
normal <- binned.plot(felm.est = normal.fe,
             noci = TRUE,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "Mean Monthly Max. Temp. in \u2103",
             ylabel = "Percentage Point Change in Pr(Physically Active)",
             panel = "a",
             panel.x = -4.5,
             panel.y = -11.5,
             linecolor = 'gray40',
             pointfill = 'gray40')

## overweight
ov <- binned.plot(felm.est = ov.fe,
             noci = TRUE,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "",
             y.axis.theme = element_blank(),
             y.axis.line = element_blank(),
             y.axis.ticks = element_blank(),
             panel = "",
             panel.x = -4.5,
             panel.y = -11.5,
             linecolor = 'gray70',
             pointfill = 'gray70')
weight <- normal + ov$layers[[2]] + ov$layers[[3]]

## obese
ob <- binned.plot(felm.est = ob.fe,
             noci = TRUE,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "",
             y.axis.theme = element_blank(),
             y.axis.line = element_blank(),
             y.axis.ticks = element_blank(),
             panel = "",
             panel.x = -4.5,
             panel.y = -11.5,
             linecolor = 'gold',
             pointfill = 'gold')
weight <- weight + ob$layers[[2]] + ob$layers[[3]]

## create legend for plot from dummy data
dum <- data.frame(point = sort(rep(letters[1:3], 10)), x = runif(30), y = runif(30))
weight <- weight + geom_point(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
            scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("gray40", "gray70","gold"), 
                           values=c("gray40", "gray70", "gold"),
                           labels=c("Normal BMI", "Overweight", "Obese")) +
            theme(legend.title=element_blank(),
                  legend.position = c(.7, 0.2),
                  legend.text=element_text(size=16, face="bold"),
                  legend.key = element_rect(colour = "black", fill="white"))
  
## weight plot
weight <- ggMarginal(weight,
           data = brfss, 
           x = brfss[, avg.tmax], 
           type = "histogram",
           margins = "x",
           size = 20,
           xparams = list(bins = 30, fill = "gold", color = 'gray60', alpha=0.75, size=0.1))

## young
young <- binned.plot(felm.est = young.fe,
             noci = TRUE,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "Mean Monthly Max. Temp. in \u2103",
             ylabel = "Percentage Point Change in Pr(Physically Active)",
             panel = "b",
             panel.x = -4.5,
             panel.y = -11.5,
             linecolor = 'gray40',
             pointfill = 'gray40')

## medium
medium <- binned.plot(felm.est = middle.fe,
             noci = TRUE,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "",
             y.axis.theme = element_blank(),
             y.axis.line = element_blank(),
             y.axis.ticks = element_blank(),
             panel = "",
             panel.x = -4.5,
             panel.y = -11.5,
             linecolor = 'gray70',
             pointfill = 'gray70')
age <- young + medium$layers[[2]] + medium$layers[[3]]

## old
old <- binned.plot(felm.est = old.fe,
             noci = TRUE,
             plotvar = "cuttmax",
             breaks = 1,
             omit = c(28,29),
             ylimit = c(-11.5,1.5),
             xlabel = "",
             y.axis.theme = element_blank(),
             y.axis.line = element_blank(),
             y.axis.ticks = element_blank(),
             panel = "",
             panel.x = -4.5,
             panel.y = -11.5,
             linecolor = 'darkorchid4',
             pointfill = 'darkorchid4')
age <- age + old$layers[[2]] + old$layers[[3]]

## create legend for plot from dummy data
dum <- data.frame(point = sort(rep(letters[1:3], 10)), x = runif(30), y = runif(30))
age <- age + geom_point(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
            scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("gray40", "gray70","darkorchid4"), 
                           values=c("gray40", "gray70", "darkorchid4"),
                           labels=c("Under 40", "40 to 65","65 and Over")) +
            theme(legend.title=element_blank(),
                  legend.position = c(.7, 0.2),
                  legend.text=element_text(size=16, face="bold"),
                  legend.key = element_rect(colour = "black", fill="white"))
  
## age plot
age <- ggMarginal(age,
           data = brfss, 
           x = brfss[, avg.tmax], 
           type = "histogram",
           margins = "x",
           size = 20,
           xparams = list(bins = 30, fill = "gold", color = 'gray60', alpha=0.75, size=0.1))

grid.arrange(weight, age, ncol = 2)
```

```{r figure3, echo=TRUE, fig.width=22, fig.height=7, dev='png', dpi=1200, fig.cap="**Main text Figure 3**"}
# splined forecast plot

## knots for linear splines
for (i in seq(from=-4.5, to=40.5, by=5)) {
  brfss <- brfss[, paste0("avg.tmax",i) := avg.tmax - (i)]
  brfss <- brfss[avg.tmax < i, paste0("avg.tmax",i) := 0]
}
setnames(brfss, names(brfss), gsub("-", "neg", names(brfss)))

## use demeanlist from lfe package to take out fixed effects to be able to use predict.lm
demean <- demeanlist(mtx = brfss[, .SD, .SDcols=(names(brfss) %like% "exerany2|avg.tmax") & !(names(brfss) %like% "prism")],
                     fl = list(cutprcp = brfss[, cutprcp],
                               cutcloud = brfss[, cutcloud],
                               cuthumid = brfss[, cuthumid],
                               cuttrange = brfss[, cuttrange],
                               mmsa=brfss[, as.factor(mmsa)],
                               time=brfss[, as.factor(time)],
                               mmsa.season=brfss[, as.factor(paste0(mmsa,season))]))

## conduct the regression on the demeaned data
cov <- grep("avg.tmax", names(demean), value=TRUE)
cov <- paste0(cov, collapse=" + ")
fmla <- as.formula(paste("exerany2", "~", cov))
history <- lm(fmla, data=demean)

## knots for linear splines on bcsd
for (i in seq(from=-4.5, to=40.5, by=5)) {
  bcsd <- bcsd[, paste0("avg.tmax",i) := avg.tmax - (i)]
  bcsd <- bcsd[avg.tmax < i, paste0("avg.tmax",i) := 0]
}
setnames(bcsd, names(bcsd), gsub("-", "neg", names(bcsd)))

## use avg.tmax splines from bcsd data to fit values
forecast <- data.table(predict(history, bcsd[, .SD, .SDcols=names(bcsd) %like% "avg.tmax"], interval="confidence"))
bcsd$forecast <- forecast[, fit]
bcsd$lwr <- forecast[, lwr]
bcsd$upr <- forecast[, upr]

## scale by 1,000 individuals (multiply by 10 since exerany is 0/100)
bcsd <- bcsd[, ':=' (forecast = forecast*10,
             lwr = lwr*10,
             upr = upr*10)]

## clean up bcsd
for (i in names(bcsd)[!(names(bcsd)=="avg.tmax") & names(bcsd) %like% "avg.tmax"]) {
  bcsd <- bcsd[, (i) := NULL]
}

## need to remove years from model names
bcsd <- bcsd[, model := gsub("2010", "", model)]
bcsd <- bcsd[, model := gsub("2050", "", model)]
bcsd <- bcsd[, model := gsub("2099", "", model)]

## create wide version of forecast
bcsd <- bcsd[, year := year(time)] # create year variable
bcsd <- bcsd[, month := month(time)] # create month variable
bcsd <- bcsd[, doy := as.numeric(strftime(time, format = "%j"))] # create doy term for casting

## cast data to wide format by year
bcsd <- dcast.data.table(bcsd, mmsa + model + month + doy ~ year, 
                         value.var = c("tmax","avg.tmax","forecast","lwr","upr"))

## get differences due to climate change between 2010 and 2050 and 2050 and 2099
bcsd <- bcsd[, diff.2050 := forecast_2050-forecast_2010]
bcsd <- bcsd[, diff.2050.lwr := lwr_2050-lwr_2010]
bcsd <- bcsd[, diff.2050.upr := upr_2050-upr_2010]
bcsd <- bcsd[, diff.2099 := forecast_2099-forecast_2010]
bcsd <- bcsd[, diff.2099.lwr := lwr_2099-lwr_2010]
bcsd <- bcsd[, diff.2099.upr := upr_2099-upr_2010]

## get mean for each month
diff.plot <- bcsd[, .(diff.2050 = mean(diff.2050, na.rm=T),
                      diff.2050.lwr = mean(diff.2050.lwr, na.rm=T),
                      diff.2050.upr = mean(diff.2050.upr, na.rm=T),
                      diff.2099 = mean(diff.2099, na.rm=T),
                      diff.2099.lwr = mean(diff.2099.lwr, na.rm=T),
                      diff.2099.upr = mean(diff.2099.upr, na.rm=T)
                      ),
                        by=.(model,month)]

mean.plot <- diff.plot[, .(mean.2050 = mean(diff.2050, na.rm=T),
                      mean.2099 = mean(diff.2099, na.rm=T),
                      hi2050 = max(diff.2050.upr, na.rm=T),
                      lo2050 = min(diff.2050.lwr, na.rm=T),
                      hi2099 = max(diff.2099.upr, na.rm=T),
                      lo2099 = min(diff.2099.lwr, na.rm=T)),
                  by=.(month)]

## plot mean and min/max (linerange) of forecasts across all models for each month
## 2050
plotforecast2050 <- ggplot(data=mean.plot) +
    geom_bar(aes(x=month, y=mean.2050, width=0.8), color="gray70", fill="gold", stat="identity") +
    geom_linerange(aes(ymax = hi2050, ymin = lo2050, x=month), color="gray40") +
    coord_cartesian(ylim = c(-12,20), xlim = c(0.5,12.5)) +
    scale_x_continuous(breaks=c(1:12), labels=c('J','F','M','A','M','J','J','A','S','O','N','D')) +
    ylab("Predicted Additional Physically Active People\n(per 1,000 Individuals)") +
    xlab("Months in 2050") +
    annotate("text", label="c", x = 0.5, y = -12, size=15) +
    geom_hline(yintercept = 0) +
    theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
          plot.title = element_text(hjust=0.025),
          axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
          axis.title.y = element_text(size=20, angle=90, vjust=1, face="bold"),
          axis.text.y = element_text(size=18, angle=0, vjust=.5, face="bold"),
          axis.text.x = element_text(size=18, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(),
          legend.position="none"
    )

## 2099
plotforecast2099 <- ggplot(data=mean.plot) +
    geom_bar(aes(x=month, y=mean.2099, width=0.8), color="gray70", fill="darkorchid4", stat="identity") +
    geom_linerange(aes(ymax = hi2099, ymin = lo2099, x=month), color="gray40") +
    coord_cartesian(ylim = c(-12,20), xlim = c(0.5,12.5)) +
    scale_x_continuous(breaks=c(1:12), labels=c('J','F','M','A','M','J','J','A','S','O','N','D')) +
    ylab("Predicted Additional Physically Active People\n(per 1,000 Individuals)") +
    xlab("Months in 2099") +
    annotate("text", label="d", x = 0.5, y = -12, size=15) +
    geom_hline(yintercept = 0) +
    theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
          plot.title = element_text(hjust=0.025),
          axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
          axis.title.y = element_text(size=20, angle=90, vjust=1, face="bold"),
          axis.text.y = element_text(size=18, angle=0, vjust=.5, face="bold"),
          axis.text.x = element_text(size=18, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(),
          legend.position="none"
    )

## plot of temperature in 2010, 2050, and 2099

## get averages across all mmsa
temp.plot <- bcsd[, .(temp.2010 = mean(tmax_2010, na.rm=T),
                      temp.2050 = mean(tmax_2050, na.rm=T),
                      temp.2099 = mean(tmax_2099, na.rm=T)),
                  by=.(model, month)]

mean.temp.plot <- bcsd[, .(temp.2010 = mean(tmax_2010, na.rm=T),
                      temp.2050 = mean(tmax_2050, na.rm=T),
                      temp.2099 = mean(tmax_2099, na.rm=T)),
                  by=.(month)]

tempforecast <- ggplot() +
                  geom_line(data=mean.temp.plot, aes(x=month, y = temp.2010, color="gray40"), alpha=1, linetype='solid', size=1.5) +
                  geom_line(data=mean.temp.plot, aes(x=month, y = temp.2050, color="gold"), alpha=1, linetype='solid', size=1.5) +
                  geom_line(data=mean.temp.plot, aes(x=month, y = temp.2099, color="darkorchid4"), alpha=1, linetype='solid', size=1.5) +
                  annotate("text", label="b", x = 1, y = -2, size=15) +
                  coord_cartesian(ylim = c(-2,40), xlim = c(1,12)) +
                  scale_x_continuous(breaks=c(1:12), labels=c('J','F','M','A','M','J','J','A','S','O','N','D')) +
                  ylab("Mean Max. Temp. in \u2103") +
                  xlab("Months") +
                  scale_color_manual(limits=c("gray40", "gold", "darkorchid4"), values=c("gray40", "gold", "darkorchid4"),labels=c("2010", "2050", "2099")) +
                  guides(fill = guide_legend(override.aes = list(colour = NULL))) +
                  theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
                        plot.title = element_text(hjust=0.025),
                        axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
                        axis.title.y = element_text(size=20,angle=90, vjust=1, face="bold"),
                        axis.text.y = element_text(size=18,angle=0, vjust=.5, face="bold") ,
                        axis.text.x = element_text(size=18, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line = element_line(),
                        legend.title=element_blank(),
                        legend.position = c(.58, 0.25),
                        legend.text=element_text(size=16, face="bold"),
                        legend.key = element_rect(colour = "black")
                  )

## historical tmax histogram
tmax.hist <- ggplot(data = brfss, aes(x = avg.tmax)) +
                  geom_histogram(aes(y=..density..), fill='gray70', color='gray40', alpha=1) +
                  ylab("Density") +
                  xlab("Mean Max. Temp. in \u2103") +
                  geom_vline(xintercept = 28.5, color='darkorchid4', linetype=2, size=1.25) +
                  coord_cartesian(xlim=c(-11.25, 46.25),ylim = c(-0.002,.0425)) +
                  scale_x_continuous(breaks=c(-10,0,10,20,30,40)) +
                  annotate("text", label="a", x = -11.25, y = -0.002, size=15) +
                  theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
                        plot.title = element_text(hjust=0.025),
                        axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
                        axis.title.y = element_text(size=20,angle=90, vjust=1, face="bold"),
                        axis.text.y = element_text(size=18,angle=0, vjust=.5, face="bold") ,
                        axis.text.x = element_text(size=18, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line = element_line(),
                        legend.title=element_blank()
                  )

grid.arrange(tmax.hist, tempforecast, plotforecast2050, plotforecast2099, ncol = 4)
```

```{r figure4, echo=TRUE, fig.width=8, fig.height=4, dev='png', dpi=1200, fig.cap="**Main text Figure 4**"}
# map forecast annual

## get state shapefiles
setwd("./shp/")
us <- invisible(readOGR(dsn=getwd(), layer="cb_2014_us_state_5m")) #read in each spatial polygon data frame
us <- us[!(us$NAME=="Puerto Rico" | us$NAME=="United States Virgin Islands" | us$NAME=="American Samoa" | us$NAME=="Guam" | us$NAME=="Commonwealth of the Northern Mariana Islands" | us$NAME=="Alaska" | us$NAME=="Hawaii"),] #continental us
projection(us) <- projection(sumdiff2050)

## multiply raster values by 10 to get 1,000 individuals for forecast values
## exerany2 is 0/100
forecast2050 <- sumdiff2050*10
forecast2099 <- sumdiff2099*10

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## stack these together
forecast <- stack(forecast2050, forecast2099)
id <- c('2050','2099') # set year id
forecast <- setZ(forecast, id) # set year as z variable
names(forecast) <- c('2050','2099') # set names

## plot raster forecasts
theme = rasterTheme(region=c('red4', 'red3','red', 'gold', 'goldenrod3', 'gray60','gray40', 'darkorchid2','darkorchid3','darkorchid4', 'purple4')) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(2,1), # horizontal layout
               names.attr=c("2050","2099"), # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Predicted Annual Change in Physically Active Person-Months Due to Climate Change\n(per 1,000 individuals)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(us, lwd=0.8, col="black")) # add us shapefiles

p
```

```{r figure5, echo=TRUE, fig.width=8, fig.height=8, dev='png', dpi=1200, fig.cap="**Main text Figure 5**"}
# map forecast by month

## get state shapefiles
setwd("./shp/")
us <- invisible(readOGR(dsn=getwd(), layer="cb_2014_us_state_5m", verbose=FALSE)) #read in each spatial polygon data frame
us <- us[!(us$NAME=="Puerto Rico" | us$NAME=="United States Virgin Islands" | us$NAME=="American Samoa" | us$NAME=="Guam" | us$NAME=="Commonwealth of the Northern Mariana Islands" | us$NAME=="Alaska" | us$NAME=="Hawaii"),] # continental us
projection(us) <- projection(sumdiff2050)

## multiply raster values by 10 to get 1,000 individuals for forecast values
## exerany2 is 0/100
forecast2050.1 <- diff2050.1*10
forecast2050.2 <- diff2050.2*10
forecast2050.3 <- diff2050.3*10
forecast2050.4 <- diff2050.4*10
forecast2050.5 <- diff2050.5*10
forecast2050.6 <- diff2050.6*10
forecast2050.7 <- diff2050.7*10
forecast2050.8 <- diff2050.8*10
forecast2050.9 <- diff2050.9*10
forecast2050.10 <- diff2050.10*10
forecast2050.11 <- diff2050.11*10
forecast2050.12 <- diff2050.12*10

forecast2099.1 <- diff2099.1*10
forecast2099.2 <- diff2099.2*10
forecast2099.3 <- diff2099.3*10
forecast2099.4 <- diff2099.4*10
forecast2099.5 <- diff2099.5*10
forecast2099.6 <- diff2099.6*10
forecast2099.7 <- diff2099.7*10
forecast2099.8 <- diff2099.8*10
forecast2099.9 <- diff2099.9*10
forecast2099.10 <- diff2099.10*10
forecast2099.11 <- diff2099.11*10
forecast2099.12 <- diff2099.12*10

## stack these together
forecast2050 <- stack(forecast2050.1, forecast2050.2, forecast2050.3, forecast2050.4, forecast2050.5, forecast2050.6, forecast2050.7, forecast2050.8, forecast2050.9, forecast2050.10, forecast2050.11, forecast2050.12)

forecast2099 <- stack(forecast2099.1, forecast2099.2, forecast2099.3, forecast2099.4, forecast2099.5, forecast2099.6, forecast2099.7, forecast2099.8, forecast2099.9, forecast2099.10, forecast2099.11, forecast2099.12)

id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(forecast2050, id2050) # set year as z variable
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(forecast2099, id2099) # set year as z variable
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
theme = rasterTheme(region=c('red4', 'red3', 'red2', 'gold3', 'gold2', 'gold1', 'gray70','gray50', 'darkorchid1', 'darkorchid3', 'darkorchid4', 'purple4')) # define color palette

nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Predicted Change in Physically Active Persons Due to Climate Change\n(per 1,000 individuals)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(us, lwd=0.8, col="black")) # add us shapefiles

p
```
